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Abstract 

We  formulate  a  Multi-Element  generalized  Polynomial  Chaos  (ME-gPC)  method  to 
deal  with  long-term  integration  and  discontinuities  in  stochastic  differential  equa¬ 
tions.  We  first  present  this  method  for  Legendre-chaos  corresponding  to  uniform 
random  inputs,  and  subsequently  we  generalize  it  to  other  random  inputs.  The 
main  idea  of  ME-gPC  is  to  decompose  the  space  of  random  inputs  when  the  rela¬ 
tive  error  in  variance  becomes  greater  than  a  threshold  value.  In  each  subdomain 
or  random  element,  we  then  employ  a  generalized  Polynomial  Chaos  expansion.  We 
develop  a  criterion  to  perform  such  a  decomposition  adaptively,  and  demonstrate 
its  effectiveness  for  ODEs,  including  the  Kraichnan-Orszag  three-mode  problem,  as 
well  as  advection-diffusion  problems.  The  new  method  is  similar  to  spectral  element 
method  for  deterministic  problems  but  with  h-p  discretization  of  the  random  space. 
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1  Introduction 


Polynomial  chaos  is  a  non-statistical  approach  to  represent  randomness  and 
is  based  on  the  homogeneous  chaos  theory  of  Wiener  [1] .  In  its  original  form 
a  spectral  expansion  was  employed  based  on  the  Hermite  orthogonal  polyno¬ 
mials  in  terms  of  Gaussian  random  variables.  This  expansion  was  applied  by 
Ghanem  and  co-workers  to  various  problems  in  mechanics  [2,3].  A  broader 
framework,  called  “generalized  Polynomial  Ghaos  (gPG)”,  was  introduced  in 

*  Corresponding  author. 

Email  address:  xlwan@dain.brown.edu  and  gk@dain.brown.edu  (Xiaoliang  Wan 
and  George  Em  Karniadakis). 


Preprint  submitted  to  Elsevier  Science 


9  March  2005 


Report  Documentation  Page 

Form  Approved 

0MB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 

1.  REPORT  DATE 

09  MAR  2005 

3.  DATES  COVERED 

00-03-2005  to  00-03-2005 

4.  TITLE  AND  SUBTITLE 

An  Adaptive  Multi-Element  Generalized  Polynomial  Chaos  Method  for 
Stochastic  Differential  Equations 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Brown  University, Division  of  Applied  Mathematics, 182  George 

Street, Providence, RI, 02912 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distrihution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE 

unclassified  unclassified  unclassified 

31 

Standard  Form  298  (Rev.  8-98} 

Prescribed  by  ANSI  Std  Z39-18 


[4,5].  This  extension  includes  a  family  of  orthogonal  polynomials  (the  so-called 
Askey  scheme)  from  which  the  trial  basis  is  selected,  and  can  represent  non- 
Gaussian  processes  more  efficiently;  it  includes  the  classical  Hermite  polyno¬ 
mial  chaos  as  a  subset.  For  example,  uniform  distributions  are  represented  by 
Legendre  polynomial  functionals,  exponential  distributions  by  Laguerre  poly¬ 
nomial  functionals,  etc.  The  method  includes  also  discrete  distributions  with 
corresponding  discrete  eigenfunctions  as  trial  basis;  e.g.,  Poisson  distributions 
are  represented  by  Charlier  polynomial  functionals. 

More  specihcally,  stochastic  ordinary  differential  equations  (ODEs)  were  con¬ 
sidered  in  [4]  and  gPC  was  shown  to  exhibit  exponential  convergence  in  ap¬ 
proximating  stochastic  solutions  at  hnite  (early)  times.  However,  the  absolute 
error  may  increase  gradually  in  time  and  become  unacceptably  large  for  long¬ 
term  integration.  Increasing  the  polynomial  order  adaptively  can  somewhat 
alleviate  this  problem,  however,  the  stochastic  solution  may  become  increas¬ 
ingly  complicated,  which  may  give  rise  to  serious  computational  difficulties. 
For  example,  if  the  stochastic  solutions  are  periodic  with  random  frequencies, 
gPC  will  lose  its  effectiveness  rapidly  due  to  the  amplihed  phase  shift  with 
time.  The  same  is  true  for  time-dependent  simulations  of  fluid  flows,  which 
are  the  problems  considered  in  [5].  In  addition,  for  discontinuous  dependence 
of  the  solution  on  the  input  random  data,  gPC  may  converge  slowly  or  fail  to 
converge  even  in  short-time  integration.  This  situation  represents  essentially  a 
discontinuity  of  the  approximated  solution  in  random  space,  for  which  global 
solutions  converge  slowly.  Therefore,  more  efficient  and  robust  schemes  are 
needed  to  enhance  the  performance  of  generalized  as  well  as  the  original  poly¬ 
nomial  chaos.  To  this  end,  a  new  method,  termed  the  Wiener-Haar  method, 
was  developed  in  [6,7]  based  on  wavelets;  its  primary  aim  was  to  address  prob¬ 
lems  related  to  the  aforementioned  discontinuities  in  random  space. 

In  this  paper,  we  develop  a  simple  but  effective  scheme  based  on  gPC,  i.e.,  we 
maintain  a  spectral  polynomial  trial  basis.  It  is  motivated  by  two  observations: 

(1)  gPC  is  more  efficient  for  relatively  small  degree  of  random  perturbation, 
and 

(2)  most  of  the  statistics  we  are  interested  in,  such  as  mean  and  variance,  are 
dehned  as  integrations  involving  the  probability  density  function  (PDF). 

To  this  end,  we  decompose  the  space  of  random  inputs  into  small  elements. 
Subsequently,  in  each  element  we  generate  a  new  random  variable  and  apply 
gPC  again.  Since  the  degree  of  perturbation  in  each  element  is  reduced  pro¬ 
portionally  to  the  size  of  random  elements,  we  can  maintain  a  relative  low 
polynomial  order  for  gPC  in  each  element.  This  multi-element  gPC  method 
(ME-gPC)  can  achieve  h-p  convergence  (as  in  spectral  elements  for  spatial 
discretization),  where  h  is  determined  by  the  size  of  random  elements  and  p  is 
the  polynomial  chaos  order.  The  concept  of  h-convergence  used  in  this  work  is 
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similar  to  that  in  [8],  where  the  basis  of  the  standard  hnite  element  method  is 
employed.  In  ME-gPC,  orthogonal  basis  (Legendre-chaos)  is  used  in  each  ran¬ 
dom  element  for  efficiency.  By  extension,  we  can  say  that  in  [6,7]  the  concept  of 
h-convergence  is  also  used  with  h  representing  the  number  of  resolution  levels 
of  wavelets.  From  the  implementation  standpoint,  the  simplicity  of  ME-gPC 
is  particularly  attractive;  for  example,  we  do  not  have  to  change  the  existing 
gPC  solver  except  for  a  subroutine  for  the  decomposition  of  random  space. 
As  we  shall  see  in  this  paper,  however,  the  results  are  dramatically  improved 
compared  to  global  gPC  expansions. 

This  paper  is  organized  in  the  following  way.  In  the  next  section,  we  recall 
the  basic  concepts  and  properties  of  gPC.  Then,  we  introduce  the  ME-gPC 
algorithm  and  the  criterion  of  the  decomposition  of  random  space  in  section  3. 
In  section  4,  we  study  the  properties  of  ME-gPC  numerically  for  several  typical 
ODE  and  PDE  problems,  including  the  open  Kraichnan-Orszag  three-mode 
problem.  A  summary  is  included  in  section  5. 


2  Generalized  Polynomial  Chaos 


The  original  polynomial  chaos  formulation  was  proposed  by  N.  Wiener  [1]. 
It  employs  Hermite  polynomials  in  terms  of  Gaussian  random  variables  as 
the  trial  basis  to  represent  stochastic  processes.  According  to  the  theorem  of 
Cameron  and  Martin  [9]  such  expansions  converge  for  any  second-order  pro¬ 
cesses  in  the  L2  sense.  The  gPC  extension  was  proposed  in  [5]  and  employs 
more  types  of  orthogonal  polynomials  from  the  Askey  scheme.  It  is  a  general¬ 
ization  of  the  Wiener’s  Hermite-chaos  and  can  deal  with  non-Gaussian  random 
inputs  more  efficiently. 

Let  (D,  JF,  P)  be  a  complete  probability  space,  where  D  is  the  sample  space, 
T  is  the  (T-algebra  of  subsets  of  D,  and  P  is  a  probability  measure.  A  general 
second-order  random  process  X[iS)  G  T2(D,  P,  P)  can  be  expressed  by  gPC  as 

00 

=  (1) 

i=0 

where  co  is  the  random  event  and  $1(^(0;))  are  polynomial  functionals  of  degree 
p  in  terms  of  the  multi-dimensional  random  variable  ^  =  (,^1, . . . ,  The 
family  is  an  orthogonal  basis  in  L2{^1,  P,  P)  with  orthogonality  relation 

=  (4>2)d,,,  (2) 

where  Sij  is  the  Kronecker  delta,  and  (•,  •)  denote  the  ensemble  average.  Here, 
the  ensemble  average  can  be  dehned  as  the  inner  product  in  the  Hilbert  space 
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in  terms  of  the  random  vector  ^ 


{/(«).9(«))  =  fni)3(i)w{m 

(3) 

imigm  = 

(4) 

4 


in  the  discrete  case,  where  w(^)  denotes  the  weight  function. 


For  a  certain  random  vector  the  gPC  basis  can  be  selected  in  such  a 
way  that  its  weight  function  has  the  same  form  as  the  probability  distribution 
function  of  The  corresponding  type  of  polynomials  and  their  associated 
random  variable  ^  can  be  found  in  [4], 


3  Multi-Element  generalized  Polynomial  Chaos  (ME-gPC) 


In  this  section,  we  develop  the  scheme  of  ME-gPC  to  maintain  the  high  accu¬ 
racy  of  gPC  for  long-term  integration  and  to  resolve  effectively  discontinuities 
in  random  space. 


3. 1  Decomposition  of  Random  Space 


Let  ^  =  (,^i(a;),  .^2(1^),  •  ■  ■  denote  a  d- dimensional  random 

vector  dehned  on  the  probability  space  (hi,  JF,  P),  where  fi  are  identical  inde¬ 
pendent  distributed  (HD)  random  variables.  Here  we  assume  that  are  also 
uniform  random  variables  dehned  as  :  Q  ^  [—1, 1]  with  a  constant  PDF 

/*  =  i 

Let  D  be  a  decomposition  of  B  with  N  non-overlapping  elements 


D 


B,  =  [a^,,b^,)x[al,bl)x---x[aib% 


Pfci  =  0)  if  h  7^  k2, 


(5) 


where  k,  /ci,  /c2  =  1, 2,  •  ■ 
random  element  as 


N.  We  dehne  an  indicator  random  variable  for  each 


Zk 


1  if  ^  G  Pfc, 
0  otherwise. 


(6) 


It  is  easy  to  see  that  D  =  ^(1)  and  ^(1)  fl  z^  ^(1)  =  0  when  i  7^  j. 

Thus,  is  a  decomposition  of  the  sample  space  D.  Then,  in  each 


4 


random  element  we  define  the  following  local  random  vector  as 


C  =(Cf,  Cl,- 

subject  to  a  conditional  PDF 

1 


= 


2^FT{zk  =  1) 


where 


=  1)  =  n 


(7) 


(8) 


(9) 


2=1 


Note  that  Pr(2:fc  =  1)  >  0.  Subsequently,  we  map  C^  to  a  new  random  vector 
defined  in  [—1,1]'^ 


y = a(c‘)  =  :  %‘(i) « 1-1. 1]' 


(10) 


with  a  constant  PDF  =  {\Yi  where 


^ 


(11) 


To  this  end,  we  present  a  decomposition  of  the  random  space  of  Given  a 
system  of  differential  equations  with  random  inputs  the  output  u{^)  is  also 
measurable  on  the  probability  space  (D,  JF,  P).  Thus,  we  can  express  u{^)  in 
each  random  element  using  subject  to  a  conditional  PDF,  which  implies 
that  we  can  hrst  approximate  u{^)  locally  by  C^  on  the  probability  space 
(2:^^(1),P  n  P(-|2:^^(l))),  then  combine  all  the  information  from  each 

random  element  to  get  u{^)  in  the  whole  random  space.  Since  most  of  the 
statistics  are  integrations  with  respect  to  the  PDF,  we  do  not  have  to  guarantee 
the  absolute  continuity  in  terms  of  ^  between  random  elements.  In  other  words, 
the  following  restriction 

■*^51(0  = '*^^2(0,  C^7?inP2,  (12) 

where  Pi  and  P2  indicate  the  closure  of  two  adjacent  random  elements,  re¬ 
spectively,  is  not  required  as  in  the  deterministic  problems  since  the  measure 
of  the  interface  is  zero.  Thus,  in  random  element  k  we  can  use  gPC  locally 
to  solve  the  system  of  differential  equations  with  random  inputs  instead  of 
According  to  the  theorem  of  Cameron  and  Martin  [9],  gPC  will  converge 
to  m(C^)  in  the  L2  sense.  Hence,  we  decompose  the  original  problem  to  N 
independent  problems  in  N  random  elements. 

In  practice  we  implement  gPC  according  to  instead  of  to  take  ad¬ 
vantage  of  the  Legendre-chaos.  After  we  obtain  the  approximation 
k  =  1,  2,  •  ■  ■  ,  A^,  of  a  random  held,  we  can  reconstruct  the  m-th  moment  of 
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u{^)  on  the  entire  random  domain  by  the  Bayes’  theorem  and  the  law  of  total 
probability  [10] 


r  1  ^  r  1 

Am  ('“(«))  =  /  «  y  Pr(^i  =  1)  /  (13) 

Jb  2  ^[-1,1]''  2 

Since  we  consider  second-order  processes  in  this  work,  m  =  1,2.  For  conve¬ 
nience,  we  use  Jk  to  denote  ^T{zk  =  1)  in  the  presentation  below. 


3.2  Accuracy 


Theorem  1  Suppose  ^  is  a  random  vector  defined  on  [—1, 1]'’*  with  IID  uni¬ 
form  components.  If  the  random  space  of  ^  is  decomposed  into  N  disjoint 
elements  with  each  element  k  described  by  a  new  uniform  random  vector 
(see  equation  (10)),  the  m-th  (m  =  1,2)  moment  of  random  field  E 
L2(I2,  T ,  P)  can  be  approximated  by  Uk{^^),  k  =  1,2,  ■  ■  ■  ,N,  with  a  L2  error 

^  =  (14) 

where  is  the  local  L2  error  of  the  m-th  moment  in  random  element  k, 
Jk  =  l^r{zk  =  1)  and  Uk{^^)  is  obtained  from  gPC. 


PROOF.  Let  u{^)  be  the  approximate  random  held.  We  hrst  assume  that 
the  m-th  moment  of  u{^)  takes  the  form 

N 

fi”(«  =  E«r(ft(«)^/i-.  (15) 

k=l 

since  B  =  ^  -B*  and  ^  E  B  (see  equation  (7)  and  (11)).  Then, 


N 

= y  pr($  e  BO  /  b"‘(c'-')  -  fir(j/=«‘)))  eaE 

k=i 

^  r  2  1 

=  1:  Pr($  £  BO  /  „  b’”(90‘(?'=))  -  fi?(y ))  (;;)"d«‘ 

k=i  “'i-i-h  ^ 

N 

=  E4^^- 

k=l 
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For  the  second  step,  we  employ  the  Bayes’  theorem  and  the  law  of  total  prob¬ 
ability  [10].  If  gPC  is  employed  to  approximate  ek  goes  to  zero  ac¬ 

cording  to  the  theorem  given  by  Cameron  and  Martin  [9].  Since  Jk  =  1, 
e  also  goes  to  zero.  Note  here  that  although  we  approximate  the  random  field 
locally  we  can  rebuild  the  global  random  held  by  equation  (15).  □ 


Note  that  Y.k=i  Jk  =  1-  Thus  is  a  weighted  mean  of  e|,  /c  =  1,  2,  •  •  •  ,  N. 
From  the  transform  (11)  we  can  see  that  the  degree  of  random  perturbation 

for  each  dimension  of  ^  is  scaled  down  from  0(1)  to  0(  "  ^  " ).  This  means 
that  the  decomposition  of  random  space  can  effectively  decrease  the  degree  of 
randomness.  Thus,  for  the  same  polynomial  order  any  would  be  smaller  than 
the  error  given  by  gPC  on  the  entire  random  space  without  the  decomposition 
of  random  space. 

In  [8]  error  estimates  were  derived  for  the  mean  and  the  variance  for  a  sim¬ 
ilar  decomposition  of  random  space  in  the  framework  of  deterministic  hnite 
element  method,  as  follows 

|h  -  h|  <  la'  -  d'l  <  C'2(p)0(h'(^’+^)),  (16) 

where  the  element  size  h  oc  N~^  in  our  case  and  p  is  the  polynomial  order.  In 
[8],  the  same  basis  as  the  deterministic  hnite  element  is  employed  to  approxi¬ 
mate  the  random  held,  where  the  accuracy  mainly  relies  on  the  decomposition 
of  random  space.  In  ME-gPC,  we  employ  Legendre- chaos  locally  to  take  ad¬ 
vantage  of  orthogonality  and  related  efficiencies. 

Let  us  now  return  to  the  two  specihc  problem  we  aim  to  address  in  this  paper: 
discontinuity  and  long-term  integration.  If  a  discontinuity  exists  in  random 
space,  then  gPC  may  converge  very  slowly  or  give  rise  to  0(1)  error.  However, 
ME-gPC  can  overcome  this  difficulty.  Let  us  assume  that  the  discontinuity 
occurs  in  the  random  element  k.  From  equation  (14)  we  can  see  that  the  error 
contribution  of  element  k  is  which  is  determined  by  the  local  approx¬ 

imation  error  in  element  k  and  the  factor  together.  So  the  error  contribution 
can  be  decreased  by  the  factor  (dictated  by  the  element  size)  even  if  the 
local  approximation  error  is  big.  Thus,  we  can  maintain  a  high  accuracy  on  the 
entire  random  domain  by  using  bigger  elements  for  the  smooth  part  (p-type 
convergence)  and  smaller  elements  for  the  discontinuous  part  (h-type  conver¬ 
gence).  To  control  the  error  in  long-term  integration  problems,  one  choice  is 
to  increase  the  polynomial  order  adaptively.  However,  the  stochastic  system 
will  become  bigger,  which  may  lead  to  a  complicated  system  of  deterministic 
differential  equations  with  all  stochastic  modes  coupled  together,  especially 
in  problems  with  high-order  nonlinearity.  In  ME-gPC,  we  can  use  a  relative 
low  polynomial  order  in  each  random  element  since  the  local  degree  of  per¬ 
turbation  has  been  scaled  down;  thus,  the  complexity  is  effectively  controlled. 
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In  practice,  the  polynomial  order  cannot  be  increased  arbitrarily  high,  which 
means  that  the  range  of  application  of  gPC  is  indeed  limited.  It  is  obvious 
that  such  a  range  can  be  effectively  extended  by  the  decomposition  of  random 
space. 


3.3  Adaptive  Criterion 


Let  us  assume  that  the  gPC  expansion  of  random  held  in  element  k  is 

Np 

=  (17) 

i=0 

where  p  is  the  highest  order  of  polynomial  chaos  and  Np  denotes  the  total 
number  of  basis  modes  given  by 


{p  +  d)\ 
p\d\ 


The  approximate  global  mean  can  be  expressed  as 


N 

^  ^  '  'dkftJk- 

k=l 


(18) 


(19) 


From  the  orthogonality  of  gPC  we  can  obtain  the  local  variance  approximated 
by  polynomial  chaos  with  order  p 


i=l 


(20) 


and  the  approximate  global  variance 


Let  7fc  be  the  error  of  the 
variance  as 


We  dehne  the  local  decay 
each  element  as  follows 


=  E  [^l,P  +  (“Lo  -  uf\  Jk-  (21) 

k=l 

term  a^  p  +  {ukp  —  uY-  We  obtain  the  exact  global 


N 


a 


=  +  X)  jkJk,  ■ 


(22) 


k=l 


rate  of  relative  error  of  the  gPC  approximation  in 


Vk 


Eiv.+i '“«<*?) 


(J. 


k,p 


(23) 


For  /i-type  refinement,  we  consider  two  factors:  the  decay  rate  of  the  relative 
error  rjk  in  each  element  and  the  factor  Jt-  We  will  split  a  random  element 
into  two  equal  parts  when  the  following  condition  is  satisfied 


vtJk  >0i,  0  <  a  <  1, 


(24) 


where  a  is  a  prescribed  constant. 

When  the  random  elements  become  smaller,  (i.e.,  Jk  becomes  smaller),  the 
value  of  rjk  satisfying  the  criterion  will  be  bigger.  Thus,  the  criterion  relaxes 
the  restriction  on  the  accuracy  of  the  local  variance  for  smaller  elements  since 
the  error  contribution  of  small  random  elements  will  be  dictated  by  their  size. 
From  equation  (22)  we  can  see  that  to  achieve  a  certain  level  of  accuracy, 
say  /5,  we  need  J2k=i'lkJk/<^‘^  ~  0{(3).  However,  it  is  difficult  to  estimate 
such  a  global  error  since  it  is  related  to  both  h-type  convergence  and  p-type 
convergence.  By  noting  the  hierarchical  structure  of  orthogonal  polynomial 
chaos  basis,  we  replace  7fc/cT^  with  r]k  and  use  rjkJk  as  an  indicator  of  the  error 
contribution  of  each  element  in  this  work. 

There  are  two  reasons  to  use  the  power  of  rjk  with  respect  to  a  in  the  criterion: 

(1)  The  decomposition  of  random  space  would  terminate  when  ~  0i.  From 
the  criterion,  we  can  see  that  rjk  must  satisfy  rjk  >  {Oi/Jk)^^°'  to  trigger 
the  decomposition  of  random  space.  If  Jk  <  9i,  rjk  must  be  greater  than 
1  and  increase  quickly  as  Jk  becomes  smaller  further  by  noting  that  both 
9i/  Jk  and  1/a  are  greater  than  1.  It  is,  in  general,  hard  to  reach  such  a 
large  rjk  in  practice,  even  for  problems  involving  stochastic  discontinuities. 
Thus,  9i  acts  as  a  limit  of  the  size  of  random  elements.  In  this  paper,  we 
usually  set  a  to  be  1/2. 

(2)  In  stochastic  discontinuity  problems  the  largest  error  contribution  is  rjkJk  ~ 
0{Jk)  ~  0{9i)  because  the  relative  error  rjk  could  be  almost  0(1)  in  the 
elements  containing  discontinuities.  For  such  a  case,  we  have  to  keep  the 
error  contribution  of  0(^1)  because  it  is  the  best  that  gPC  can  do;  how¬ 
ever,  we  can  eliminate  the  error  contribution  of  random  elements  with¬ 
out  discontinuities.  Note  that  rjkJk  ~  0{ril~°‘9i),  where  9i  is  weighted  by 
T]k~^-  Thus,  in  random  elements  without  discontinuities  the  error  contri¬ 
bution  will  be  much  smaller  than  9i  since  r]k  <  1  in  these  elements.  Fi¬ 
nally,  the  total  error  contribution  Y.k=i  VkJk  would  be  0{mJk)  ~  0{m9i), 
where  m  is  the  number  of  random  elements  with  0(6^1)  error  contribu¬ 
tion.  So,  rj]r^  works  as  a  filter  and  6*1  also  acts  as  an  accuracy  threshold 
besides  the  aforementioned  limit  of  element  size. 

Furthermore,  we  use  another  threshold  parameter  ^2  to  choose  the  most  sen¬ 
sitive  random  dimension.  We  dehne  the  sensitivity  of  each  random  dimension 
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as 


ri  = 


i  =  1,2,  -  ■  ■  ,d 


(25) 


where  we  drop  the  subscript  k  for  clarity  and  the  subscript  -i^p  denotes  the 
mode  consisting  only  of  random  dimension  with  polynomial  order  p.  All 
random  dimensions  which  satisfy 


r*  >  6*2  ■  max  rj,  0  <  6*2  <  1,  i  =  1,2,  ■  ■  ■  ,  d  (26) 

will  be  split  into  two  equal  random  elements  in  the  next  time  step  while 
all  other  random  dimensions  will  remain  unchanged.  Hence,  we  can  reduce 
the  total  element  number  while  gaining  efficiency.  Considering  that  h-type 
rehnement  is  efficient  in  practice,  we  only  present  results  given  by  h-type 
rehnement  in  this  work.  For  some  cases,  say  stochastic  discontinuity  problems, 
h-type  rehnement  may  be  the  most  effective  choice  since  p-type  convergence 
may  not  be  maintained  anymore.  This  is,  of  course,  not  surprising  given  what 
we  know  for  deterministic  problems  [11]. 


3.4  Numerical  Implementation 


When  h-type  rehnement  is  needed,  we  have  to  map  the  random  held  from  one 
mesh  of  elements  to  a  new  mesh  of  elements.  Suppose  that  the  gPC  expansion 
of  the  current  random  held  is 


Np 

i=0 


(27) 


then  we  assume  that  the  gPC  expansion  in  the  next  level  takes  the  following 
form 

=  w  (^  (^))  (28) 

i=0 

where  ^  G  [—1,  l]'^.  To  determine  the  {Np  +  1)  coefficients  Ui,  we  choose  {Np  +  1) 
points  i  =  0, 1,  •  •  ■  ,  Np,  which  are  the  uniform  grid  points  in  [—1, 1]^  and 
solve  the  following  linear  system 


"hoo  "hio  ■  ■  ■  ‘hvpO 

Uo 

$01  $11  ■  ■  ■  ^Npl 

Ui 

= 

(j-‘  («,)) 

^ONp  ^INp  ■  ■  ■  ^NpNp 

(29) 
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where  <l>jj  =  <l>j  .  We  rewrite  the  above  equation  in  matrix  form  as 

Au  =  u.  (30) 

Due  to  the  hierarchical  structure  of  the  basis,  exists  for  any  {Np  +  1) 
distinct  points  in  [—1, 1]'^.  When  h-type  rehnement  is  implemented  we  divide 
the  random  space  of  a  certain  random  dimension  into  two  equal  parts.  For 
example,  if  corresponds  to  element  [a,  b]  in  the  original  random  space  [—1, 1], 
the  elements  [a,  and  [^,b]  will  be  generated  in  the  next  level.  However, 
due  to  the  linearity  of  transformation,  we  do  not  have  to  perform  such  a  map 
from  the  original  random  space,  as  we  can  just  separate  the  random  space  of 
which  is  [—1, 1],  to  [—1,0]  and  [0, 1].  Therefore,  the  matrix  A  will  be  the 
same  for  every  h-type  rehnement,  and  we  only  need  to  compute  A~^  once  and 
store  it  for  future  use.  When  rehnement  is  needed,  we  can  obtain  u  easily  by 
a  matrix-vector  multiplication 


u  =  A  ^u.  (31) 

For  a  relatively  small  polynomial  order  [p  <  10),  the  mapping  cost  is  small. 
Now  we  summarize  the  ME-gPC  algorithm. 

Algorithm  1 

-Step  1:  construct  a  stochastic  ODE/PDE  system  by  gPC 
-Step  2:  perform  the  decomposition  of  random  space  adaptively 
time  step  i:  from  1  to  A 

loop  over  all  random  elements 
ifp'^Jk  >  di,  then 

if  f'n  >  ^2  ■  then 

split  random  dimension  fn  into  two  equal  ones  and  generate 
local  random  variables  fn,i  and  fn,2 
end  if 
end  if 

map  information  to  the  children  random  elements 
update  the  information  of  new  elements  by  gPC 
end  loop 
end  time  step 

-Step  3:  postprocessing  stage 


3.5  Generalization 


Let  C  be  a  general  (i.e.,  non-uniform)  random  vector,  whose  components  are 
IID  random  variables.  Let  (  denote  any  component  of  C.  We  can  approximate 
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it  by  Legendre-chaos  in  the  form 


C  =  (32) 

i=0 

where  is  a  uniform  random  variable.  The  procedure  for  such  an  approxi¬ 
mation  can  be  found  in  [4],  Note  here  that  we  need  d  IID  uniform  random 
variables  to  approximate  all  components  of  C-  By  expressing  everything  in 
terms  of  the  Legendre-chaos,  then  we  can  employ  ME-gPC  in  terms  of 


Another  choice  is  to  first  decompose  the  random  space  of  (.  Assume  that  u{() 
is  a  random  held  of  then  the  m-th  moment  of  m(C)  is 

=  /  u^iOHOdC,  (33) 

Jb 

where  /i(C)  is  the  PDF.  Suppose  that  we  have  decomposed  the  random  space 
of  C  to  elements  5*,  i  =  1,  2,  •  ■  ■  ,  iV.  The  above  equation  can  be  rewritten  as 

f^m{u)=J2[  ViU^iO^^dC,  (34) 

JBi  Vi 

where  Vi  =  /g,  h{C)dC.  We  can  then  express  as  a  conditional  PDF  of  (  in 

B, 

h{C\Bi)  =  (35) 

Vi 

Then,  the  m-th  moment  of  m(C)  can  be  expressed  in  the  following  form 


N 


=  U^iOHCmdC. 


2=1 


(36) 


Now  we  can  employ  the  hrst  choice  to  approximate  the  conditional  PDF 
h{C,\Bi)  by  uniform  random  variables  f.  Since  we  approximate  h{C,\Bi)  only 
in  a  subspace  of  C,  we  may  use  a  smaller  number  of  Legendre-chaos  modes  for 
a  desired  level  of  accuracy. 


Finally,  another  choice  is  to  construct  orthogonal  polynomials  on-the-hy  for 
arbitrary  PDFs.  This  construction  is  under  development  (see  [12]). 


4  Numerical  Results 


In  this  section,  we  hrst  demonstrate  the  convergence  of  MF-gPC  for  an  al¬ 
gebraic  equation  and  a  simple  ODF.  Next,  we  focus  on  issues  related  to  dis¬ 
continuities  in  random  space  and  study  the  Kraichnan-Orszag  problem.  Sub¬ 
sequently,  we  present  numerical  results  for  the  stochastic  advection-dihusion 
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equation.  Finally,  we  demonstrate  the  h-type  convergence  of  the  decomposi¬ 
tion  of  random  space  for  the  approximation  of  general  random  inputs. 


4-1  A  Simple  Algebraic  Equation 


We  hrst  revisit  the  following  stochastic  algebraic  equation  considered  in  [8] 

cu  =  1,  (37) 

where  c  is  a  positive  uniform  random  variable  in  [a,  b]. 


In  Fig.  1,  the  h-type  convergence  is  shown,  with  the  mean  on  the  left  and 
variance  on  the  right.  Here  we  set  a  =  2  and  6  =  3.  By  a  least-squares  £t  of 
the  data,  we  obtain  that  the  index  of  algebraic  convergence  is  2(p-|- 1)  for  both 
the  mean  and  the  variance,  which  is  consistent  with  the  theoretical  estimates 
given  in  [8]. 


Fig.  1.  /i-type  convergence  for  the  algebraic  equation.  Left:  Mean;  Right:  Variance. 


4-2  One- Dimensional  ODE 


In  this  section  we  study  the  performance  of  ME-gPC  for  the  following  simple 
ODE  equation  studied  with  the  original  gPC  in  [4] 

d-ti 

—  =  —k{uj)u,  m(0;ci;)  =  uo,  (38) 

where  k,{u)  ~  U{—1, 1).  The  exact  solution  can  be  easily  found  as 

(39) 

In  Fig.  2,  we  show  the  exponential  convergence  of  ME-gPC  for  different 
meshes.  We  can  see  that  for  greater  number  of  equidistant  random  elements. 
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not  only  is  the  error  smaller,  but  the  rate  of  convergence  is  much  sharper.  We 
show  the  algebraic  convergence  of  ME-gPC  in  terms  of  element  number  N  in 
Fig.  3.  For  this  problem,  the  algebraic  index  of  convergence  is  2(p  +  1)  for 
both  mean  and  variance,  which  means  e  O  (^N  2(p+i)^_  yve  have  obtained  a 
large  algebraic  index  of  convergence,  which  implies  that  random  elements  can 
influence  the  accuracy  dramatically.  In  Fig.  4,  the  error  evolution  of  gPC  and 
MF-gPC  is  shown  for  two  different  levels  of  accuracy.  Because  the  accuracy  of 
exact  solutions  is  set  to  be  10“^'^,  there  is  some  oscillation  at  the  beginning  of 
the  curves.  It  can  be  seen  that  when  the  error  of  gPC  becomes  big  enough,  6i 
can  trigger  the  decomposition  of  random  space  and  the  accuracy  can  then  be 
improved  signihcantly.  In  Fig.  5,  we  show  how  the  number  of  random  elements 
increases  adaptively.  Note  here  that  the  mesh  can  be  non-uniform,  because  we 
only  decompose  the  random  elements  in  which  the  criterion  is  satisfied. 


Fig.  2.  Stochastic  ODE:  Exponential  convergence  of  ME-gPC  with  respect  to  poly¬ 
nomial  order  (t  =  5).  Left:  Mean;  Right:  Variance. 


Fig.  3.  Stochastic  ODE:  Algebraic  convergence  of  ME-gPC  with  respect  to  number 
of  random  elements  {t  =  5).  Left:  Mean;  Right:  Variance. 
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Fig.  4.  Stochastic  ODE:  Error  evolution  of  gPC  and  ME-gPC  with  a  =  1/2.  Left: 
Mean;  Right:  Variance. 


Fig.  5.  Stochastic  ODE:  Error  (left  axis)  and  number  of  random  elements  (right 
axis)  with  a  =  1/2.  Left:  p  =  3,6i  =  0.01;  Right:  p  =  3,0i  =  0.001. 


4-3  The  Kraichnan-Orszag  Three-Mode  Problem 


It  is  well  known  that  polynomial  chaos  fails  in  a  short  time  for  the  so-called 
Kraichnan-Orszag  three- mode  problem  [13].  In  this  section  we  hrst  explain 
why  this  happens  and  subsequently  we  apply  ME-gPC  to  effectively  resolve 
this  fourty-year  old  open  problem. 


4-3.1  Why  gPC  fails 

The  Kraichnan-Orszag  problem  [13]  is  a  nonlinear  three-dimensional  stochas¬ 
tic  ODE  system: 
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(40) 


We  first  check  the  deterministic  solutions  of  equation  (40).  Given  different 
initial  conditions,  deterministic  solutions  can  be  basically  separated  into  four 
different  groups  Qi,  i  =  1,2,  3, 4,  which  are  shown  in  Fig  6.  All  these  four 
groups  of  solutions  are  periodic.  If  the  initial  conditions  are  located  on  the 
planes  xi  =  X2  and  xi  =  —X2,  the  corresponding  solutions  would  stay  on 
these  two  planes  forever  due  to  two  hxed  points  (0,  0,  \j2xi{Q)  +  x^O))  and 

(0,  0,  —^2x\{Q)  +  xKO)).  By  considering  the  properties  of  elliptic  functions 
[14],  we  can  obtain  the  analytic  solutions  of  each  group.  Here  we  only  give  the 
analytic  form  of  group  gi\ 

Xi  =  Pcn[g(t  -  to)],  xs  =  Qdn[g(t -to)],  Xg  = -Psn[g(t  -  to)],  (42) 

where  cn[-],  sn[-]  and  dn[-]  are  Jacobi’s  elliptic  functions  and  P,  Q,  P,  q  and  to 
are  constants  to  be  determined.  We  now  substitute  equation  (42)  into  equation 
(40)  to  obtain 

Pq  =  QR,  Qk'^q  =  PR,  Rq  =  2PQ,  (43) 

where  k  is  the  modulus  of  elliptic  functions.  Since  we  have  three  initial  condi¬ 
tions 


Pcn[g(t  -  to)]  =a;i(0;a;), 

Qdn[q{t  -  to)]=X2{0;u), 

-Rsn[q{t  -  to)]=X3{0;u),  (44) 

we  have  six  equations  with  six  unknowns  P,  Q,  R,  k,  q  and  to.  Thus,  we  have 
obtained  the  exact  general  solution  of  the  Kraichnan-Orszag  problem. 


We  now  consider  the  following  initial  conditions 

a;i(0)  =a  +  0.01^,  0:2(0)  =  1.0,  0:3(0)  =  1.0,  (45) 

where  is  a  uniform  random  variable  and  a  is  a  constant.  By  solving  equation 
(43)  and  (44),  we  can  determine  the  unknowns  as 
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(46) 


P^  =  fiO  +  l,  =  |  i?'  =  2f(0  +  l, 

p2  =  3^  /^2  ^  ‘^fiO  +  ^,  to  =  -dn"^[^]/p, 

where  f{^)  =  a  +  0.01.^. 

Next  we  examine  the  Fourier  expansions  of  Jacobi’s  functions: 


sn[u] 

cn[M] 

dn[u] 


2tt 

kK 

2% 

kK 

71 


sin  q 

- + 

1-q 

cos  2: 


,3/2 


sin  3z  q 
+  - 


sin  5z 


l  +  q 


+ 


1  —  q^ 
cos  3z 


Q 


271 

2K^1{ 


l  +  q^ 

.2, 


+ 


Q 


1  —  q^ 
cos  5z 


+  ■■■ 


1  + 
.3, 


+ 


qcos2z  q^  cos  4:Z  q'^  cos6z 

+  ^  +  - ^  +  - 


1  + 


l  +  g^ 


1  +  g® 


(47) 


where  g  =  q{^),  K  =  K{^),  and  =  z{^,  t).  First,  we  can  see  that  the  frequency 
depends  on  the  random  variable  It  is  well  known  that  this  will  reduce  the 
effectiveness  of  gPC  as  the  initial  phase  difference  will  be  amplihed  very  fast 
as  time  increases.  In  Fig.  7,  we  show  how  the  period  of  xi  change  as  a;i(0)  — 1. 
We  can  see  that  the  period  of  xi  will  increase  to  inhnity  as  a;i(0)  goes  to  1. 
Note  here  that  if  a;i(0)  =  1,  the  initial  point  (1, 1, 1)  would  be  on  the  plane 
Xi  =  X2.  Second,  if  g  goes  to  1,  it  is  clear  that  we  need  more  and  more  terms 
for  the  expansion  of  sn[u],  which  means  that  the  order  of  polynomial  chaos 
must  increase  correspondingly  to  resolve  the  solution. 


From  equations  (45)  and  (46)  we  can  see  that  if  ^  is  uniform  in  [—1,1],  xi 
is  uniform  in  [a  —  0.01,  a  +  0.01]  and  the  range  (non-uniform)  of  k{^)  is 
[^1(0;  —  0.01)  -|-  i,  +  0.01)  -|-  |].  Let  kr  denote  the  upper  bound  of  k{^). 

It  is  clear  that  if  a  — 0.99,  kr  1.  By  the  properties  of  elliptic  functions, 
we  know  that  g  — 1  when  k  ^  1.  Thus,  for  the  same  degree  of  perturbation 
gPC  should  work  less  efficiently  when  a  0.99,  because  k{^)  will  be  closer 
to  1.  Now,  we  investigate  four  simple  cases:  a  =  0.94,  0.96,  0.98,  and  0.99. 
For  simplicity  we  only  show  the  results  for  xi,  since  the  situation  is  similar 
for  X2  and  x^,-  In  Fig.  8  we  show  how  gPC  fails  when  a  0.99.  It  can  be  seen 
that  in  Fig.  8  (a)-(d)  the  valid  range  of  polynomial  chaos  with  order  p  =  6 
becomes  shorter  as  a  increases.  If  a  is  strictly  less  than  0.99  corresponding  to 
g  <  1,  increasing  the  polynomial  order  can  efficiently  improve  the  results  of 
polynomial  chaos.  For  the  cases  (a)-(c),  the  results  of  polynomial  chaos  with 
order  p  =  20  agree  very  well  with  the  results  of  Monte  Carlo  with  100,  000  re¬ 
alizations.  However,  if  a  =  0.99  is  included,  the  periods  of  stochastic  solutions 
will  change  from  a  hnite  value  to  inhnity  and  increasing  the  polynomial  order 
hardly  improves  the  results  for  this  case.  It  is  shown  in  (d)  that  the  correct 
part  of  the  variance  given  by  polynomial  chaos  with  order  p  =  30  is  almost 
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the  same  with  that  given  by  polynomial  chaos  with  order  p  =  6.  Therefore,  it 
is  at  the  bifurcation  point  where  gPC  fails  to  converge. 

In  general,  if  the  initial  random  data  does  not  intersect  with  the  planes  xi  =  X2 
and  Xi  =  —X2,  we  can  improve  the  results  of  polynomial  chaos  by  increasing 
the  polynomial  order,  otherwise,  polynomial  chaos  will  diverge  even  after  a 
short  time  of  integration. 


Fig.  6.  Deterministic  solutions  of  the  Kraichnan-Orszag  problem  subject  to  different 
initial  conditions.  Left;  3D  phase  space;  Right;  2D  projection  on  X1-X2  plane. 


■  0  5  10  15  20  25  30  35  40  45  50 


Fig.  7.  Kraichnan-Orszag  problem;  Several  deterministic  solutions  of  xi  versus  time 
corresponding  to  different  initial  conditions. 


4-3.2  One- Dimensional  Random  Input 

Let  us  first  study  the  random  discontinuity  of  the  Kraichnan-Orszag  three¬ 
mode  problem,  which  is  introduced  by  one-dimensional  random  input.  For 
computational  convenience  and  clarity  in  the  presentation  we  hrst  perform 
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Fig.  8.  Comparison  of  variance  obtained  from  gPC  and  Monte  Carlo  simulations, 
(a):  a  =  0.94;  (b):  a  =  0.96;  (c):  a  =  0.98;  (d);  a  =  0.99. 

the  following  transformation 


As  a  result,  we  will  rotate  the  deterministic  solutions  by  |  around  to  X3  axis 
in  the  phase  space.  Now  the  new  system  is 


dy2 

dt 


-y2ys 


dys 

dt 


-yl  +  yl 


(49) 
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subject  to  initial  conditions 

l/i(0)  =  |/i(0;o^),  1/2(0)  =  1/2(0;  0^),  1/3(0)  =  1/3(0;  0^).  (50) 

From  now  on,  we  will  study  this  problem  based  on  equation  (49).  Note  that  the 
discontinuity  occurs  at  the  planes  //i  =  0  and  //2  =  0  after  the  transformation. 
Gaussian  random  variables  are  used  as  random  inputs  in  [13].  Here,  we  use 
uniform  random  variables  since  the  discontinuity  can  be  introduced  similarly. 
Thus,  we  study  the  stochastic  response  subject  to  the  following  random  input 

yi{0;u)  =  l,  1/2(0;  cn)  =  O.l^(cn),  //3(0;cn)  =  0,  (51) 

where  f  ~  U{—1, 1).  Since  the  random  initial  data  y2{0]uj)  can  cross  the  plane 
y2  =  0,  we  know  from  the  aftermentioned  discussion  that  gPC  will  fail  for  this 
case. 

In  Fig.  9,  we  show  the  evolution  of  the  variance  of  yi  within  the  time  interval 
[0,30].  For  comparison  we  include  the  results  given  by  gPC  with  polynomial 
order  p  =  30.  It  can  be  seen  that  comparing  to  the  results  given  by  Monte 
Carlo  with  1,  000,  000  realizations,  gPC  with  polynomial  order  p  =  30  begins 
to  lose  accuracy  at  t  ~  8  and  fails  beyond  this  point  while  ME-gPC  converges 
as  61  decreases.  In  Table  1,  we  show  the  maximum  normalized  error  of  the 
variance  of  yi,  y2  and  ps  at  t  =  30  given  by  ME-gPC  and  the  corresponding 
number  of  random  elements.  It  is  seen  that  when  the  threshold  parameter  9i 
decreases,  the  accuracy  becomes  better  and  we  can  obtain  almost  0{9i)  error. 
As  we  mentioned  before,  the  reason  that  errors  are  usually  bigger  than  9i 
is  due  to  the  discontinuity  which  can  reduce  the  convergence  of  gPC.  It  can 
be  seen  that  for  the  same  polynomial  order  we  need  more  random  elements 
to  get  a  better  accuracy;  on  the  other  hand,  for  the  same  9i  increasing  the 
polynomial  order  can  reduce  the  number  of  random  elements. 

In  Fig.  10,  we  show  four  adaptive  meshes.  We  can  see  that  around  the  point 
.^  =  0  in  random  space  of  where  the  discontinuity  occurs,  the  random 
elements  are  smallest,  which  means  that  the  discontinuity  can  be  captured  by 
small  random  elements.  In  Fig.  11,  we  show  the  errors  of  Monte  Carlo  and 
ME-gPC  in  terms  of  computational  cost.  The  error  is  the  Loo  error  of  the 
variance  of  yi  in  the  time  interval  [8,  30],  where  gPC  fails.  To  implement  gPC, 
we  need  to  apply  Galerkin  projection  onto  the  chaos  basis,  resulting  in  the 
ensemble  average  of  three  basis  modes.  Here,  we  count  the  operations 

of  for  ME-gPC  in  order  to  estimate  its  cost.  For  Monte  Carlo,  the 

number  of  realizations  is  employed  in  the  cost  evaluation.  Let  n  denote  the 
number  of  operations.  If  the  data  in  Fig.  11  are  approximated  by  a  hrst-order 
polynomial  in  a  least-squares  sense,  we  can  obtain  accuracy  proportional  to 
^-0.49,  .^^-2.25^  ^-2.99  respectively,  for  Monte  Carlo  and  ME-gPC 

with  polynomial  order  p  =  3,  p  =  4  and  p  =  5,  respectively.  The  decay  rate 
for  Monte  Carlo  is  about  as  expected.  Comparing  to  Monte  Carlo,  the 
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errors  of  ME-gPC  show  a  much  greater  decay  rate  in  terms  of  the  cost.  We  can 
see  that  the  speed-up  increases  for  higher  accuracy,  which  implies  that  ME- 
gPC  is  an  efficient  alternative  to  Monte  Carlo  for  integration  where  high-order 
accuracy  is  required.  In  Fig.  12,  we  show  the  error  contribution  of  each  random 
element.  Here  we  compare  two  criteria  with  a  =  |  and  a  =  |.  It  is  seen  that 
the  shape  of  error  distribution  is  like  an  isosceles  triangle,  i.e.,  a  “Gibbs-like” 
behavior.  On  the  apex  of  the  triangle  is  the  largest  error  contribution,  where 
discontinuity  occurs.  The  error  contribution  decreases  quickly  away  from  the 
discontinuity,  since  ~  hfc  is  much  smaller  on  the  smooth  part. 

Because  gPC  loses  accuracy  as  time  increases,  the  error  contribution  of  each 
element  will  become  larger  with  time  and  more  random  elements  with  relative 
errors  of  0(1)  would  appear  around  the  discontinuity  point.  For  a  smaller  a, 
the  error  contribution  near  the  discontinuity  decreases  much  faster. 


Fig.  9.  Evolution  of  the  variance  of  yi  for  one-dimensional  random  input. 


Table  1 

Maximum  normalized  errors  of  the  variance  of  yi,  y2  and  y^  at  t  =  30  with  a  =  1/2. 
(The  results  given  by  ME-gPC  with  6i  =  10“^  and  polynomial  order  p  =  5  are  used 
as  exact  solutions.). 


01  =  10-2 

01  =  10-3 

01  =  10-^ 

01  =  10-5 

N 

Error 

N 

Error 

N 

Error 

N 

Error 

p  =  3 

46 

3.10e-2 

106 

2.32e-3 

280 

1.37e-4 

820 

2.87e-5 

p  =  4 

36 

9.90e-2 

74 

3.24e-3 

138 

3.45e-4 

286 

2.31e-5 

p  =  5 

28 

7.24e-2 

44 

4.10e-3 

78 

2.90e-4 

130 

4.35e-6 

21 


(b) 


(d) 


Fig.  10.  Adaptive  meshes  for  the  ID  random  input  with  a  =  1/2.  (a);  Q\  =  0.01, 
p  =  3;  (b):  6»i  =  0.001,  p  =  3;  (c):  6»i  =  0.0001,  p  =  3;  (d):  6»i  =  0.0001,  p  =  5. 

J^-.3.3  Two-Dimensional  Random  Input 

In  this  section  we  use  ME-gPC  to  study  the  Kraichnan-Orszag  problem  with 
two-dimensional  random  input 

yi{0;u)  =  1,  P2(0;cu)  =  O.l^i(cu),  i/3(0;cu)  =  6(cu),  (52) 

where  and  .^2  are  uniform  random  variables  with  unit  standard  deviation. 

In  Fig.  13,  we  show  the  evolution  of  the  variance  of  yi,  p2  and  y^  and  an 
adaptive  two-dimensional  mesh.  For  comparison  we  include  the  result  given  by 
gPC  with  polynomial  order  p  =  10.  It  can  be  seen  that  gPC  with  polynomial 
order  p  =  10  begins  to  diverge  around  t  ~  4  while  MF-gPC  with  p  =  5 
Legendre-chaos  shows  good  convergence  to  the  results  given  by  Monte  Carlo 
with  1,  000,  000  realizations.  From  the  final  refined  mesh,  we  can  see  that  the 
results  are  more  sensitive  to  .^1,  because  .^1  can  cross  the  plane  P2  =  0  where  the 
discontinuity  occurs.  Note  here  that  the  discontinuity  domain  is  a  line.  In  Fig. 
14,  we  show  the  error  of  Monte  Carlo  and  MF-gPC  in  terms  of  computational 
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Fig.  11.  Error  verus  cost  of  Monte  Carlo  simulations  and  ME-gPC  with  different 
polynomial  orders  (based  on  the  L^o  error  of  the  variance  of  yi  in  the  time  interval 
[8,  30]).  Here  we  only  count  the  average  number  of  operations  in  one  time  step. 

cost.  Here  we  regard  the  results  given  by  ME-gPC  with  9i  =  10“®  and  p  =  5 
as  exact  solutions.  From  the  empirical  fit  we  obtain  an  accuracy  proportional 
to  and  respectively,  for  Monte  Carlo  and  ME-gPC  with 

p  =  3  and  p  =  5.  It  is  seen  that  ME-gPC  is  much  faster  than  Monte  Carlo 
for  higher  accuracy.  Comparing  to  the  ID  case,  however,  the  decay  rate  of 
relative  error  becomes  smaller  because  both  the  random  dimension  and  the 
discontinuity  domain  become  larger. 


4-3.4  Three-Dimensional  Random  Input 

In  this  section  we  use  ME-gPC  to  study  the  Kraichnan-Orszag  problem  with 
three-dimensional  random  input 

l/i(0)  =  6(^),  1/2(0)  =  6(^),  1/3(0)  =  6(^),  (53) 

where  .^1,  .^2  and  .^3  are  uniform  random  variables  with  unit  standard  deviation. 

In  Fig.  15,  we  show  the  evolution  of  variance.  Due  to  the  symmetry  of  yi  and 
1/2  in  equation  (49)  and  the  symmetry  of  |/i(0)  and  1/2(0)  in  the  random  inputs, 
the  variances  of  yi  and  y2  are  the  same.  Here  we  only  show  the  results  for  yi 
and  ys-  It  can  be  seen  that  gPC  diverges  around  t  ~  1  and  fails  subsequently 
while  ME-gPC  shows  good  convergence  as  before.  For  this  case,  the  random 
space  [—1,1]^  of  random  inputs  contains  both  yi  =  0  and  y2  =  0  where 
discontinuities  occur.  Comparing  to  the  case  with  2D  random  inputs,  the 
discontinuity  domain  is  much  larger.  Thus,  it  is  more  difficult  to  resolve  the 
3D  case.  Based  on  the  results  given  by  ME-gPC  with  polynomial  order  3  and 
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(b) 


Index 


Index 


(d) 


Fig.  12.  Error  contribution  of  each  random  element  given  by  two  criteria  with  dif¬ 
ferent  a.  01  =  10“^  and  p  =  5.  (a):  a  =  1/2, t  =  50;  (b);  a  =  1/2, t  =  100;  (c): 
a  =  1/4,  t  =  50;  (d):  a  =  1/4,  t  =  100. 


di  =  10“^,  the  Loo  errors  of  the  variance  of  yi  in  the  time  interval  [1.5,  6]  are 
0.16%  and  0.21%,  respectively,  for  Monte  Carlo  with  1,000,000  realizations 
and  ME-gPC  with  polynomial  order  p  =  3  and  6i  =  10“^.  Thus,  these  two 
errors  are  comparable.  For  this  case,  the  speed-up  of  ME-gPC  is  much  lower 
compared  to  the  2D  problem.  From  the  previous  results,  we  know  that  this 
speed-up  would  increase  for  higher  accuracy,  but  the  increasing  speed  would  be 
lower  comparing  to  the  ID  and  2D  cases.  In  Fig.  16,  we  show  the  evolution  of 
the  random  elements  generated.  It  can  be  seen  that  to  maintain  the  accuracy, 
the  element  number  has  to  increase  at  a  speed  about  100  elements  per  time 
unit. 

In  summary,  ME-gPC  shows  good  convergence  when  solving  the  Kraichnan- 
Orszag  problem  and  it  can  achieve  a  desired  accuracy  at  a  cost  much  lower 
than  Monte  Carlo.  However,  ME-gPC  loses  efficiency  for  problems  with  strong 
discontinuity  and  high-dimensional  random  inputs,  because  the  number  of 
random  elements  has  to  increase  fast  to  maintain  a  desired  accuracy. 
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Fig.  13.  The  Kraichnan-Orszag  problem  with  2D  random  inputs,  a  =  1/2, 
01  =  0.1,0.01,0.001  and  62  =  0.1.  (a);  a^;  (b):  a^;  (c):  ■  (d):  adaptive  mesh  for 

01  =  0.001  and  p  =  5. 

4-4  Stochastic  Advection- Diffusion  Equation 


In  this  section  we  consider  the  2D  stochastic  advection-diffusion  equation  first 
studied  in  [15]  using  gPC 

^(x,  +  u(x;a;)  ■  V0  =  z/V^0  (54) 

where  u(x;a;)  =  {y  +  a{u),  —x  —  b{u)).  For  the  initial  condition 

0(x,  0;  00)  =  e-[(— o)^+fe-yo)^]/2A2^  ^55^ 

the  corresponding  exact  solution  can  be  found  as 
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Fig.  14.  Error  versus  cost  of  Monte  Carlo  simulations  and  ME-gPC  with  2D  Legen- 
dre-Chaos  (based  on  the  Loo  error  of  the  variance  of  yi  in  the  time  interval  [4, 10]). 
Here  we  only  count  the  average  number  of  operations  in  one  time  step. 


Eig.  15.  Evolution  of  variance  for  the  3D  Kraichnan-Orszag  problem.  02  =  10  ^ 
Left:  Right: 

where  A  is  a  constant  and 


{X  =  X  +  b{uj)  —  {xq  +  b{uj))  cost  —  {yo  +  a{uj))  sint, 
y  =  y  +  a{uj)  +  {xq  +  b{uj))  sint  —  (yo  +  a{uj))  cost. 

Here  we  let  a{u)  =  b{u)  =  0.1.^,  where  ^  ~  t/(— 1,1).  In  Fig.  17,  we  show 
the  convergence  of  ME-gPC  with  eqnidistant  elements,  p-type  convergence  on 
the  left  and  h-type  convergence  on  the  right.  We  can  see  that  ME-gPC  not 
only  exhibits  exponential  converge  but  shows  an  increasing  convergence  rate 
as  the  number  of  elements  increases.  For  h-type  convergence,  we  only  show 
the  results  of  up  to  four  random  elements,  since  the  error  decreases  quickly. 
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Fig.  16.  Evolution  of  the  element  number  for  the  3D  Kraichnan-Orszag  problem. 
01  =  10-3. 

It  is  seen  that  the  index  of  algebraic  convergence  is  related  to  the  polynomial 
order,  where  the  decay  rate  corresponding  to  higher  polynomial  order  is  very 
large.  More  experiments  are  required  to  estimate  the  exact  convergence  rate 
numerically. 


Fig.  17.  Convergence  of  gPC  and  ME-gPC  for  2D  advection-diffusion  equation  with 
a{uj)  =  b{(jj)  =  0.1^  att  =  3.14.  Left:  p-type  convergence;  Right:  h-type  convergence. 


4-5  Approximation  of  a  Beta-type  Random  Variable  by  Legendre- Chaos 

Finally,  we  demonstrate  how  to  generalize  ME-gPC  to  other  random  variables. 
We  consider  a  Beta-type  random  variable  Y  of  distribution  ti^e{a,j3),  where 
i^e(Q;,  j3)  is  the  conventional  dehnition  of  Beta  distribution  in  the  domain  [0, 1] 
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Here  B{-,  ■)  denotes  the  Beta  function.  Let  a  =  1  and  (3  =  0,  then  the  PDF 
of  Y  is 

f{y)  =  ‘2y.  (58) 

Since  the  uniform  random  variable  used  in  Legendre-chaos  is  dehned  in  the 
domain  [—1,1],  we  introduce  a  new  random  variable  X  dehned  in  [—1, 1]  with 
the  transformation  Y  =  Thus,  the  PDF  of  X  is 

fix)  =  (59) 

Let  us  assume  that  the  random  space  [—1,1]  of  X  is  separated  into  N  equal 
elements  [a,b].  In  each  element  we  dehne  a  new  random  variable  Xi,i  = 
1,  2,  ■  ■  ■  ,N  with  a  corresponding  PDF 


fiiXi)  = 


l+Xi 


l+Xi 


I[a,b]  /(^)dt 


(1  +  a/2  +  h/2){h  —  a) 


(60) 


Subsequently,  we  use  a  uniform  random  variable  r  to  express  Xj.  A  transfor¬ 
mation  of  variables  in  probability  space  shows  that 


-dr  =  fi{xi)dxi  =  dF{xi), 


(61) 


where  F  is  the  distribution  function  of  Xj.  Thus,  we  can  obtain 


=  Fax  (62) 

After  inverting  the  above  equation,  we  obtain 

Xi  =  F~^{^  "^  '^)  =  ^(1  -F  a/2  +  b/2){b  -  a)(l  r)  (1  a)^  -  1.  (63) 

Then  Xj  can  be  expressed  by  Legendre-chaos  as 


with 


p 

Xi  =  J2xi,j^jir) 


i=i 


1  +  r 
2 


|<*,(T)-dT. 


(64) 


(65) 


Now  each  Xj  has  been  approximated  by  a  uniform  random  variable  r;  thus,  we 
can  implement  ME-gPC  in  each  element  when  solving  a  stochastic  differential 
equation  with  random  inputs  related  to  X.  Here,  we  only  check  the  accuracy 
of  /i2(X)  =  E[X^].  We  compute  /i2(X)  using  equation  (36).  In  Fig.  18,  we 
show  the  error  of  /i2(X)  in  terms  of  the  element  number  N.  It  is  seen  that  an 
algebraic  convergence  with  index  —4  is  obtained,  which  means  that  the  error 
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is  proportional  to  This  specific  value  is  dictated  by  the  accuracy  of  the 
mapping  that  we  performed  and  can  be  improved  if  higher  accuracy  is  desired. 
Therefore,  the  decomposition  of  random  space  can  also  be  used  to  approximate 
a  general  random  variable  in  order  to  improve  accuracy.  Furthermore,  we  can 
use  a  low-order  Legendre-chaos  when  implementing  ME-gPC  in  each  random 
element. 


Fig.  18.  Error  of  /U2(k)  for  a  Beta  distribution. 


5  Summary 

We  have  extended  the  generalized  polynomial  chaos  (gPC)  framework,  first 
presented  in  [4,5],  to  a  multi-element  formulation  (ME-gPC).  The  new  ap¬ 
proach  can  maintain  a  desired  accuracy  by  adaptively  decomposing  the  ran¬ 
dom  space  of  random  inputs  when  a  simple  criterion  is  satished.  Correspond¬ 
ingly,  the  efficiency  and  especially  the  effectiveness  of  gPC  is  signihcantly 
improved. 

To  investigate  the  performance  of  ME-gPC  we  present  several  examples  in¬ 
cluding  stochastic  algebraic,  ordinary  and  partial  differential  equations.  In 
particular,  we  address  errors  in  long-time  integration  and  in  discontinuities  in 
random  space.  An  example  with  one-dimensional  ODE  shows  that  ME-gPC 
can  achieve  h  —  p  type  of  convergence.  The  error  of  long-term  integration  is 
efficiently  controlled  by  the  criterion  we  developed  for  the  adaptive  decom¬ 
position  of  random  space.  Subsequently,  we  explain  why  gPC  fails  for  the 
classical  Kraichnan-Orszag  three-mode  problem,  and  study  it  with  ME-gPC 
for  different  random  inputs.  The  results  indicate  that  ME-gPC  can  capture 
accurately  the  discontinuity  by  the  decomposition  of  random  space.  In  partic¬ 
ular,  the  adaptive  criterion  can  be  used  to  select  the  most  sensitive  random 
dimension,  and  thus  make  the  decomposition  of  random  space  more  efficient. 
A  two-dimensional  advection-diffusion  equation  is  also  simulated  by  ME-gPC. 
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The  results  suggest  that  ME-gPC  could  also  improve  the  efficiency  of  gPC  for 
stochastic  PDEs.  More  results  for  stochastic  problems  of  incompressible  flow 
using  the  ME-gPC  method  presented  here  are  included  in  [12],  Finally,  we 
approximate  a  random  variable  of  Beta  distribution  by  Legendre- Chaos,  thus 
demonstrating  how  to  deal  with  general  non-uniform  random  inputs. 

ME-gPC  is  efficient  for  stochastic  systems,  which  contain  no  or  small  subdo¬ 
mains  of  discontinuities,  such  as  the  ID  ODE  model  and  the  Kraichnan-Orszag 
problem  with  ID  or  2D  random  inputs.  However,  its  efficiency  is  reduced 
signihcantly  by  the  rapidly  increasing  number  of  random  elements  for  prob¬ 
lems  with  high-dimensional  random  inputs  and  large  discontinuities,  as  in  the 
Kraichnan-Orszag  problem  with  3D  random  inputs.  Such  problems  require 
new  approaches  in  constructing  appropriate  low-dimensional  approximations, 
as  in  the  work  of  [16,17]. 
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